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ADDENDUM  TO  REPORT  ON  THE 
OIL-RPM  COMPUTER  CODE 


This  addendum  describes  a  number  of  improvements  to  the  version  of 
OIL-RPM  documented  in  GAMD-84-97.  In  EDIT  and  MAP  the  changes  represent 
primarily  a  streamlining  of  the  program.  The  small  change  to  CDT  ensures 
stability  in  a  more  rational  manner.  The  requirements  for  enlarging  the 
grid  are  also  described,  and  a  flow  chart  of  the  equation  of  state  sub¬ 
routine  (ES)  is  Included. 
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REVISED  VERSIONS  OF  CDT,  EDIT  AND  MAP 

The  improved  version  of  CDT  computes  the  time  step  using  the  speed 
|u|  +  C  of  a  sound  wave  relative  to  the  grid,  where  u  denotes  the  maxi¬ 
mum  of  the  radial  and  axial  velocities.  This  allows  the  user  to  set  the 
stability  fraction  (STAB)  to  .6,  or  perhaps  larger,  and  retains  stability; 
where  in  the  earlier  version  it  was  recommended  that  STAB  be  set  to  .4. 
This  version  of  CDT  can  be  directly  substituted  for  the  one  listed  pre¬ 
viously. 

However,  to  use  the  improved  versions  of  EDIT  and  MAP  the  following 
statements  must  be  added  to  subroutine  INPUT  after  statement  #50: 

FRSTD  -  1. 

if(kunit.le.o)  kunit  *  7 

IF(MSYMBL.LE.O)  MSYMBL  •  26 

The  prominent  features  of  the  new  EDIT  routine  are; 

(1)  It  prints  the  angular  distributions  of  the  mass,  momenta, 
and  energy. 

(2)  It  uses  a  variable  (KUNIT)  instead  of  a  constant  to  specify 
the  unit  number  of  the  dump  tape.  (The  user  can  start  a 
second  dump  tape  by  setting  KUNIT  a  9  in  the  restart  deck. 

The  program  will  read  from  unit  7‘  but  will  vrite  out  on 
unit  9.) 

(3)  An  end  of  file  mark  is  written  at  the  end  of  each  restart 
dump.  (The  next  dump  writes  over  this  end  of  file  mark.) 

(4)  The  general  flow  and  organization  of  EDIT  have  been  simplified 
making  the  routine  easier  to  modify. 


The  new  MAP  routine  is  more  efficiently  coded,  and  uses  only  one 
variable,  MSYMBL,  to  determine  the  degree  of  resolution  for  all  the  maps 


subroutine;  cot 

c  . 

c 

DIMENSION  AMX(2502) * AIX (2502) *U(2502)  *V(2502)  »P(2502)  » 

1  X(S2)  r XX (54)  »TAU(52)  *JPM(52)  * 

2  Y( 102 )  *  YY  ( 104)  »FLEFT(102)*  YAMC(102>»  SIGC(102)* 

3  GAMC ( 102) # 

4  PK ( 15) *  2(150)  * 

5  XP(26*51) * YP( 26*51 ) » 

6  PLC204)  *UL(204)  »PR(204)  * 

7  RSNC52) *  RST(52)» 

8  CMXPC5)  * CMYP (5)  *IJ(5)  *JK(5)  * 

9  DX (52)  *DDX (54)  »DY(102)  *DDY(104)  * 

S  SNB ( 52 )  *STB ( 52)  *UK(52*3)  *VK(o2*3)  *RH0(52*3) 

C  ***  DIMENSIONED  ARRAYS 

C  ***  2-BLOCK  IS  SAVED  ON  TAPE. 


COMMON 

Z 

COMMON 

PK 

COMMON 

YY* 

XX 

COMMON 

DDX* 

DDY 

COMMON 

AMX* 

AIX* 

U* 

V»  P 

COMMON 

TAU* 

JPM 

COMMON 

UL  * 

PL 

COMMON 

XP  * 

YP* 

CMXP*  CMYP 

***  NON-DIMENSIONED 

VARIABLES 

COMMON 

AID 

*  AMMV 

,ammy  ,ampy 

»amur  *amut 

* AMVR  » 

1AMVT 

»DELEB 

*DELER 

»DELET 

»delm  »DTODX 

*DXYMIN*EAMMP 

*eampy  * 

2E 

*ERDUMP* I 

*13 

»IWS  *J 

*  K  *KA 

*  KB  » 

3LL 

*MD 

*ME 

•  MZT 

»NERR  »NK 

*NPRINT» 

4NR 

*NRZ 

*NULLE 

•PIOTS 

,siemin»snr 

*SNT  »STR 

* SOLID  * 

5SUM 

*TESTRH 

*TWOPI 

*URR 

*WS  * WSA 

» WSB  »WSC 

*WFLAGF* 

6WFLAGL 

* WFLAGP 

C 

C  ***  THE  FOLLOWING  EQUIVALENCES  MAKE  AVAILABLE 

C  X(0)*  Y ( 0)  *  DX (0) »  DY(0) 

C 

EQUIVALENCE  (XX  (2) *  X(l))»  (VY(2)»  Y(D) 

EQUIVALENCE  (DDX(2)»  DX(1)),  (DDY (2)  *  DY ( 1 ) ) 

C 


C 

r 

*♦*  SPECIAL 

equivalences  for 

PH2  ONLY 

w 

EQUIVALENCE 

(ULfFLEFT)* 

(UL( 103) * YAMC) * 

r 

1 

(PL»GAMC*PR)* 

(PL( 103) fSIGC) 

C 

SPECIAL  EQUIVALENCES  FOR 

PH3  ONLY 

w 

EQUIVALENCE 

( UL#RSN) » 

1 

(PL*RST) * 

( P  *  UK ) * 

2  . 

(P( 157) »VK) * 

(P( 313) *  SNB ) » 

3 

( P ( 365) *  STB) » 

(P(417) *RH0) 

U 

c 

***  special 

equivalences  FOR 

EDIT 

C 

EQUIVALENCE  (PR(1)*  IJ)»  (PR(6)»  JK) 


C 

C  *♦*  2-STORAGE  EQUIVALENCES 

C 

EQUIVALENCE  (Z(  l)»PROB  )»(Z(  2) *CYCLE  >* 
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ooooooooooo 


1(Z(  3) »DT 

)»  IZ( 

4) 

#NUMSP 

)  * 

(Z( 

5) » 

NFRELP) » (Z ( 

6)' 

NDUMP7) 

' 

2(Z(  7) # ICSTOP 

)»(Z( 

8) 

#PIDY 

)» 

CZ( 

9) 

»TOPMU  )»(Z( 

10) 

'RTMU 

)' 

3(Z(  11 ) » STK1 

)»(Z( 

12) 

»numrez) » 

(Z( 

13) 

»ETH  )»(Z( 

14) 

'  UN14 

)' 

4 ( Z (  lb) » RHINIT 

)»(Z( 

16) 

» PROjI 

)» 

( z 

17) 

»UN17  )»(Z( 

18) 

» XMAX 

)' 

5(Z(  19) #  NZ 

)»(Zt 

20) 

»NREZ 

)  # 

(Z( 

21) 

t  AMDM  )»(Z( 

22) 

'UVMAX 

)' 

6 ( Z (  23)»UN23 

)>(Z( 

24) 

»DMlN 

)» 

(Z( 

25) 

»JSTR  )»(Z( 

26) 

»DTNA 

)' 

7(Z(  27)»CVIS 

)»(Z( 

28) 

#STK2 

) » 

(Z( 

29) 

»STEZ  )»(Z( 

30) 

'  NC 

)# 

8(Z(  31 ) » UN31 

)»(Z( 

32) 

» NRC 

)» 

(Z( 

33) 

» IMAX  )»(Z( 

34) 

'  IMAXA 

)' 

9 ( Z (  3b) » UMAX 

)»( Z( 

36) 

t  JMAxA 

)* 

CZ( 

37) 

»KMAX  )»(Z( 

38) 

'KMAXA 

) 

EQUIVALENCE 
1(Z(  39) »BOTM 

)»(Z( 

40) 

»B0TmV 

)» 

CZ( 

41) 

»NUMSPT) » (Z( 

42) 

'CZERO 

)' 

2(Z(  43)  » NUMSCA 

)»(Z( 

44) 

» PRLlM 

)* 

(Z( 

45) 

»PRDELT)» (Z( 

46) 

»PRFACT) 

EQUIVALENCE 
1(Z(  47) ' 11 

)»(Z( 

48) 

*  12 

)* 

czc 

49) 

» IPCYCL) » (Z( 

50) 

»TSTOP 

)' 

2 ( Z (  51 ) » RHOF1L 

)»(Z( 

52) 

#TARGV 

)  r 

(Z( 

53) 

»N3  )»(Z( 

54) 

» IVARDY) » 

3(Z(  55) » VT 

)»(Z( 

56) 

»N6 

)* 

(Z( 

57) 

»RTM  )»(Z( 

58) 

»RTMV 

)' 

4(Z(  59) »UN59 

)»(Z( 

60) 

» N10 

)  * 

(Z( 

61) 

»N11  )»(Z( 

62) 

'GAMMA 

)' 

5 ( Z (  63) » TOPM 

)»(Z( 

64) 

»B0TmU 

)» 

(Z( 

65) 

»SN  )»(Z( 

66) 

»topmv 

)' 

6 ( Z (  67) »PRYBOT 

)»(Z( 

68) 

♦PRYtOP) » 

(Z( 

69) 

»PRXRT  )»(Z( 

70) 

»CYCPH3) » 

7(ZC  71) »REZFCT 

)»(Z( 

72) 

i  TARGl 

)  ’ 

(Z( 

73) 

»PROJU  )»(Z( 

74) 

'BBOUND) * 

8(Z(  7b)»EVAP 

)»  (Z( 

76) 

»ECK 

)» 

(Z( 

77) 

»NECYCL) * (Z( 

78) 

'II 

) ' 

9(ZC  79) » JJ 

)»(Z< 

80) 

» NMP 

)» 

(Z( 

81) 

» Y2  )  #  ( Z  ( 

82) 

'EZPH1 

) 

equivalence 

1(Z(  83) » IVARDX 

)f  (Z( 

84) 

#  T 

)* 

(Z( 

85) 

rNMPMAX) » (Z( 

86) 

»PMIN" 

) ' 

2 ( Z (  87) » INTER 

)»(Z( 

88) 

»TAYbOT)» 

(Z( 

89) 

»TAYTOP) » (Z( 

90) 

» IEMAP 

)' 

3 ( Z (  91)»MC 

)i(Z( 

92) 

» MR 

>* 

(Z( 

93) 

'MZ  )»(Z( 

94) 

'MB 

) 

EQUIVALENCE 

1(Z(  9b)»REZ 

)»(Z( 

96) 

» NODUMP) * 

(Z( 

97) 

»UN97  )»(Z( 

98) 

'UN98 

)' 

2(Z(  99) » UN99 

)»(Z(100) 

t EVApM 

)  * 

(Z ( 101) 

»EVAPEN) » ( Z ( 102) 

»EVAPMU> » 

3(Z(1Q3) »EVAPMV 

) » ( Z ( 104 ) 

»EZPH2 

)  * 

(Z(105) 

♦SNL  ) » (Z( 106) 

»STL 

4(Z( 107) » TAXRT 

) » (Z(108) 

»IDNMAP)» 

(Z ( 109) 

» IPRMAP) » ( Z ( 110) 

'ROEPS 

)' 

5(2(111) t RHlNI 

) » ( Z  ( 1 1 2 ) 

»VINI 

)’ 

( Z ( 113) 

'FINAL  ) » ( Z ( 114) 

» IVMAP 

)' 

6(Z(115) t RHOZ 

) » ( Z  ( 1 1 6 ) 

»ESA 

)  * 

CZC117) 

»ESEZ  ) » (Z( 118) 

'ESB 

)' 

7(Z( 119) #ESCAPA 

) » (Z(120) 

*ESESP 

)* 

(Z(121) 

»ESESG  ) » (Z( 122) 

'ESES 

)' 

8(Z(123) »ESALPH 

) » ( Z ( 124 ) 

»ESBETA) * 

( Z ( 125) 

»ESC APB) ' (Z( 126) 

» IUMAP 

)' 

9(Z( 127) $  SSI 

)» (Z(128) 

»SS2 

)» 

(Z (129) 

'UMIN  ) ' (Z ( 130) 

'SS4 

) 

EQUIVALENCE 

- 

1(Z(131) tPRTIME 

)  *  ( Z ( 1 32 ) 

»EOR 

)• 

(Z( 133) 

»EOr  ) ' ( Z ( 1 34 ) 

'EOB 

)' 

2(Z(135) »EMOR 

)»<Z(136) 

»DXF 

)* 

(ZC137) 

»DYF  ) ' (Z ( 138) 

»RHOMlN) » 

3(Z( 139) t STAB) * 

(Z ( 140) t 

XIENRG) 

» 

( Z ( 141 ) 

'XKENRG) '  (Z( 142) »XTENRG) » 

4(Z(l43) »STT 

) » (Z(144) 

»DTMIN 

)» 

(Z(145) 

'TRNSFO  » (Z(146) 

'EMOT 

)' 

5(Z( 147) » JPROJ 

) » ( Z ( 148 ) 

»CNAUT 

)» 

(Z ( 149) 

'BBAR  ) ' (Z( 150) 

»EMOB 

) 

C  ♦♦♦  SPECIAL  EQUXV  FOR  ES  AND  CDT 

EQUIVALENCE  (RHOW'NULLE) 


END  OF  COMMON 


♦♦♦CHECK  COURANT  CONDITION  AND  PARTICLE  VELOCITY# 
♦♦♦RECORD  I  AND  J  OP  ZONE  WHERE  DT  IS  CONTROLLED# 
♦♦♦FIRST  CALCULATE  PRESSURES  FROM  EQ#  OF  ST# 

C 

C 


k 


JLO  TRIAL-O. 

bRATlO=10»**10 

C  ♦♦WSC  WILL  BE  MAXIMUM  U  OR  V 

WSC-G . 

DO  60  I=1»I1 
K=I+1 

DO  60  J-l» 12 

RHOW=AMX ( K ) / ( TAU ( I ) *OY ( J > ) 

CALL  ES 

C  ♦♦♦  JF  DENSITY  OF  CELL  K  IS  LESS  THAN  RH0MIN»  IT’S 

C  VELOCITY  OR  SOUNO  SPEED  IS  NOT  USED  IN  DETERMINING  DT. 

IF  (KHOW.LT.RHOMIN)  GO  TO  60 
IF  (ABS(P(K)).LT.PMIN)  P(K)=0. 

IF  (CNAUT.GT.O.)  GO  TO  20 
C 

C  ***CALCULATE  SOUND  SPEED  FOR  POLYTOPIC  GAS  WITH 

C  ♦♦♦GAMMA  EQUAL  TO  ESA+l. 

WSD=SQRT(GAMMA^ABS(P(K) )/RH0W) 

GO  TO  40 
C 

C  ♦♦♦CHECK  FOR  NEGATIVE  PRESSURE. 

20  IF  (P(K).GT.O.)  GO  TO  30 

C  ♦♦♦  NEGATIVE  PRESSURES  NOT  ALLOWED  ALONG  GRID  BOUNDARY 

C  AND  NOT  ALLOWED  ANYWHERE  UNTIL  ACTIVE  GRID  REACHES 

C  JSTR( INPUT  PARAMETER  FOR  TURNING  ON  STRENGTH 

C  CALCULATIONS). 

IF  ((IMAX.NE.1.ANO.I.EQ.IMAX)*OR.J.EQ.JMAX.OR.I2.LT.JSTR)  P(K)=0. 
C 

C  ♦♦♦PRESSURE  IS  NEGATIVE  OR  ZERO 

WSD=CNAUT 
GO  TO  40 

C 

C  ♦♦♦PRESSURE  IS  POSITIVE. 

30  WS=CNAUT+BBAR^SQRT(P(K) ) 

WSA=SQRT ( GAMMA*P ( K ) /RHOW ) 

WSD-AMAX1 ( WS » WSA ) 

C  ♦♦♦  WSB  IS  MAXIMUM  OF  RADIAL  AnD  AXIAL  VELOCITY  OF  CELL  K. 

40  WSB-AMAX1 ( ABS ( U ( K ) ) *  ABS ( V ( K ) ) ) 

C  ♦♦♦  WSC  STORES  MAXIMUM  VELOCITY  OF  CELLS  USED  TO  DETERMINE 

C  DT.  PRINTED  AS  MAXUV. 

WSC-AMAX1 ( WSC  t WSB ) 

C  ♦♦♦WSD  IS  SOUNO  SPEED  OF  CELL  K. 

WS=WSD+WSB 

IF  (WS.LE. TRIAL)  GO  TO  50 

C  ♦♦♦  TRIAL  STORES  MAXIMUM  VELOCITY  PLUS  SOULND  SPEED  OF 

C  CELL  USED  TO  DETERMINE  DT. 

TRIAL-WS 

C  ♦♦♦  CMAX  IS  SOUND  SPEED  OF  CELL  CONTROLLING  DT. 

CMAXswSO 

50  IF  (WS.LE. 0.)  GO  TO  60 

DXYMIN=AMIN1 (DX(I)tDY(J)) 

RATI0=DXYMIN/WS 

IF  (RATIO. GT.SRATIO)  GO  TO  60 

C  ♦♦♦  I  AND  J  OF  CELL  CONTROLLING  DT  STORED  IN  NlO  AND  Nil 

C  FOR  PRINTOUT. 

N10=I 

N11=J 

5 


c 

•c 

c 

69 

C 

C 

C 

65 

C 

C 

c 


70 


C 

C 

80 


C 


***  SR ATI 0  IS  SMALLEST  VALUE 
SRATlOsRATIO 


CALCULATED  FOR  RATIO. 


♦♦♦END  OF  h  J  LOOP 

KSK+IMAX 

♦♦♦  IF  TRIAL.LE.O.  THERE  IS  PROBABLY  AN  ERROR  IN  THE  INPUT 
PARAMETERS  FOR  ThE  INITIAL  VELOCITY »  ENERGY  OR  OENSITY 
OF  THE  PACKAGES. 

IF  (TRIAL.LE.O.)  GO  TO  180 


♦**  IF  FINAL. EQ.O. USE  STAB  FOR  VALUE  OF  STABILITY  FRACTION 
IF  FINAL. GT.O. USE  A  GEOMETRIC  PROGRESSION  WITH  STAB 
,  AS  THE  INITIAL  VALUE  AND  FINAL  AS  THE  FINAL  VALUE. 

IF  (FINAL. LT.O.)  GO  TO  70 
STAB»2.*STAB 


STABsAM INK STAB *FINAL) 

DTsSTAB^SRATIO 

IF  (NC.GT.O)  GO  TO  80 


IF  (UTMIN.GT.O. .OR.OTMIN.lt. 0. )  GO  TO  80 
DTMIN  =  (10.*^2) ♦DT 


♦♦♦  IS  CONTROL-CELL  ISOLATED 
K=(NU-1)^IMAX«-N10«-1 
WSSO, 

IF  (N10.GT.1)  WS=AMX(K-1) 

IF  (N10.LT.IMAX)  WSsAMX(K*l)+WS 
IF  (Nll.GT.l)  WS=AMX(K-lMAX)«-WS 
IF  (Nll.LT. JMAX)  WS=AMX(K-HMAX)*WS 
IF  (WS.GT.O.)  GO  TO  100 

♦♦♦  ISOLATEO.  SO  OESTROY  IT. 

WSs( AIX  (K)  +  (U(K) ♦♦2+V(K) ♦♦2)*»5)*AMX(K) 
EVAPM=EVAPM+AMX(K) 

EVAPENSEVAPEN+WS 

ETH=ETH-WS 

E V APMU=E VAPMU+AMX ( K ) ♦U ( K ) 

E VAPMV=E VAPMV+AMX ( K ) ♦ V ( K ) 


WRITE  (6*300)  NIOrNll 
AMX(K) SO. 

AIX(K)SO. 

P(K)S0. 

U(K)S0. 

V(K)so. 

C  ♦♦♦  RECALCULATE  DT. 

GO  TO  10 

C  ♦♦♦  INCREMENT  TIME  AND  CYCLE. 

100  T=T*DTNA 

IF  (T.LT.O.)  GO  TO  170 

NCsNC+1 

CYCLE=NC 

C  ♦♦♦  RESET  NPRINT.  NPRlNT=l  ON  PRINT  CYCLES. 

NPRINTSO 

C  ♦♦♦  OEFINE  VELOCITY  AND  ENERGY  CUTOFFS  USED  IN  MAP  AND  PH2. 

UM I N=TR I AL*ROEPS 
SIEMINSUMIN^2 
PM  I  NsRHOZ*CNAUT *UM I N 
IF  (PMIN.LT .ROEPS)  PMIN=UMIn*RHOZ^TrIaL 

WRITE  (6*310)  N10*N11*T*DT*dTMIN*CMAX*U(K) * V(K) *DX(N10) *DY(Nii ) *UM 


6 


c 

c 


104 

106 

c 

c 

c 

c 

c 

c 


110 


120 


130 


140 

150 

C 

C 

160 

C 

170 

C 

iao 

190 

c 

c 

c 

c 

c 

c 

200 


llN'PMIN 

♦  *♦  AFTER  STAb.GE. FINAL  CHECK  ON  SIZE  OF  OT.  DTMIN  CAN 
BE  DEFINED  IN  INPUT  DECK. 

IF  (STAB. LT. FINAL)  GO  TO  106 
IF  (DT.LE. DTMIN)  GO  TO  160 
CONTINUE 
DTNAsDT 

♦♦♦  TESTRH  =  ,2*RH0Z 

THE  PRESSURE  OF  COLD*FREE  SURFACE  CELLS  IS  REDUCED  BY  A 
FACTOR#F#  WHICH  ACCOUNTS  FOR  THE  EFFECT  OF  FREE  SURFACE 
LOCATION  ON  THE  PRESSURE  GRADIENT.  F  IS  THE  DENSITY  OF 
THE  LOWEST  DENSITY  ADJACENT  CELL  DIVIDED  BY  THE  NORMAL 
DENSITY*OR  F  IS  TESTRH  -  WHICHEVER  IS  SMALLEST 

wT=TESTRH 
00  150  1=1*11 
K=m 

DO  150  J=1»I2 

RHOW=AMX(K)/(DY(J)*TAU(D) 

WTB=WT 

IF  (AIX(K) . GE.ESESQ)  60  TO  150 
IF  (KHOW.LT.SOLID)  GO  TO  150 
IF  (I.EQ.IMAX)  GO  TO  110 
WTA=AMX(K+1)/(DY(J)*TAU(I+D) 

IF  (WTA.LT.WT)  WTB=WTA 
IF  (I.EQ.l)  GO  TO  120 
WTA=AMX(K-1)/(DY(J)*TAU(I-D) 

IF  (WTA.LT.WTB)  WTB=WTA 
IF  (J.EQ.JMAX)  GO  TO  130 
KA=K+IMAX 

WTA=AMX (KA)/(DY ( J+1)*TAU(I) ) 

IF  (WTA.LT.WTB)  WTB=WTA 
IF  (J.EQ.l)  GO  TO  140 
KB=K-IMAX 

WTA=AMX ( KB ) / ( DY ( J-l ) *TAU ( I ) ) 

IF  (WTA.LT.WTB)  WTB=WTA 
IF  (WTB.LT.WT)  P(K)=P(K)*WTB/RHOZ 
KSK+IMAX 
GO  TO  200 

♦♦♦  DT  TOO  SMALL 

NK=lu4 
GO  TO  190 

***  T  IS  NEGATIVE 

NK=102 
GO  TO  190 

***  DT  WILL  BE  NEGATIVE  OR  ZERO* 

NK=65 
GO  TO  190 
NR=3 

CALL  ERROR 

♦♦♦FIND  THE  MAXIMUM  PRESSURE  ON  EACH  COLUMN  AND 
♦♦♦STORE  ITS  CELL  NUMBER  AS  JPM.  THIS  WILL  BE  USED 
♦♦♦IN  DETERMINING  THE  REGION  IN  WHICH  PHASE  3  IS 
♦♦♦USED.  WSA  WILL  BE  A  RUNNING  MAXIMUM  OF  THE 
♦♦♦PRESSURE  IN  THE  GRID* 

WSAS-1.E30 
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DO  270  1=1* XI 

C  WS  WILL  BE  LOCAL  MAXIMUM  OF  COLUMN  I# 

vvS=-l.E30 
K= ( I2"l) ♦IMAX+I+1 
JP=I2 
UINTL=1 

C  ♦♦♦  START  AT  TOP  OF  COLUMN  AND  LOOK  FOR  PRESSURE  PEAK. 

210  DO  220  J=JINTL»I2 

IF  (P(K).LT.WS)  GO  TO  230 
WS=P(K) 

C  ♦*♦  JP  IS  J-XNDEX  OF  CELL  WITH  PEAK  PRESSURE. 

JP=Jp-l 

220  KSK-IMAX 

C  ♦♦♦  IF  YOU  FALL  THROUGH* THEN  THERE  WAS  NO  MAXIMUM  IN  THIS 

C  COLUMN 

GO  TO  260 
C 

C  3328  ***  COME  HERE  IF  PRESSURE  HAS  PASSED  A  LOCAL  MAXIMUM 

C  ***  PTEMP  IS  PEAK  PRESSURE  OF  COLUMN  I. 

C 

230  PTEMP=P(K+IMAX) 

IF  (PTEMP. LT.WSA)  GO  TO  240 

C  ***  WSA  WILL  BE  PEAK  PRESSURE  IN  ACTIVE  GRID  (ABSOLUTE 

C  MAXIMUM). 

WSA=PTEMP 
GO  TO  250 
C 

C  3329  ***  PTEMP  IS  LOCAL  MAXIMUM  BUT  IS  LESS  THAN  ABSOLUTE 

C  MAXIMUM 

240  IF  (PTEMP. GT.0.3*WSA)  GO  TO  250 

C 

C  ♦  ♦♦  THIS  LOCAL  MAXIMUM  IS  NOT  BIG  ENOUGH  TO  USE  FOR  UPM 

C 

JINTL=J*1 

C 

UPsUP-1  •  . 

C 

C  ***  WE  MAY  HAVE  REACHED  BOTTOM  OF  COLUMN 

IF  (UINTL.GE.I2)  GO  TO  260 

C  ***  CONTINUE  DOWN  COLUMN  SEARCHING  FOR  SUFFICIENTLY  LARGE 

C  LOCAL  MAXIMUM, 

WSSP(K) 

KSK-IMAX 
GO  TO  210 

C  **♦  IF  POSITION  OF  PEAK  PRESSURE  IN  COLUMN  I  DOES  NOT 

C  ADVANCE  FROM  ONE  CYCLE  TO  THE  NEXT*  DO  NOT  CHANGE 

C  VALUE  OF  UPM. 

250  JPsJP+1 

IF  (UP.LE.UPM(D)  GO  TO  270 
JPM(I)SJP 
C 

C  ♦  ♦♦IF  UPM  IS  ZERO  THE  SHOCK  HAS  NEVER  REACHED  THIS 

C  ♦♦♦LOCATION.  IF  IT  IS  NONZERO  THE  SHOCK  HAS  PASSED 

C  ***AND  WE  MUST  CONTINUE  TO  INCREASE  X  UNTIL  THE 

C  ♦•♦RIGHT  BOUNDARY  OF  THE  SHOCK  IS  REACHED. 

260  IF  (UPM(I).LE.O)  GO  TO  280 
C  ♦♦♦  END  OF  I  LOOP, 


e 


270  CONTINUE 

C  ***  IF  PEAK  PRESSURE  OF  COLUMN  I  HAS  GONE  BELOW  A  THIRD 

C  THE  GRID  MAXIMUM*  AND  IF  JPMUJso*  FROM  THE  PREVIOUS 

C  CYCLE*  WE  HAVE  REACHED  ThE  RIGHT  EDGE  OF  THE  SHOCK* 

280  CONTINUE 

C  ***  JPMU)  MUST  BE  MO  iOTONIC  DECREASING 

K=I1-1 

DO  290  IWSsi'K 
I=I1-IWS 

290  IF  ( JPM( I )  *LT • JPM( 1+1) )  JPM(I>=UPM<IM> 

RETURN 

C 

300  FORMAT  </4H  CDT»I3#I4»2X»3lHlS0LATED  CONTROL  CELL  DESTROYED/) 

310  FORMAT  </4H  CDT*I3#I4*4H  T=# 1PE12.6»5H  DYs* 1PE9.3»8H  DTMIN=»1PE 
19 • 3*  4H  C-* 1PE9*3»4H  U=*1PE9.3»4H  V=»1PE9.3»5H  DX=# 1PE9.3*5H  0 

2Ys»1P£9*3/13X#5HUMIN=»1PE9.3»7H  PMIN=1PE9.3> 

END 


\ 
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SUBROUTINE  EDIT 


DIMENSION  AMXC2502) #AIX<2502>  *U(2502> 


COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


X(b2)  *XX ( 54)  »TAU(52)  * 

Y (102)  *  YY ( 104)  *FLEFT(102) 

6AMC (102) * 

PK ( lb) »  Z ( 150 )  * 

XP(26»51)*YP(26»5l)» 

PL (204)  *UL(204)  *PR<204) 

RSN ( 52 ) »  RST(52)» 

C)fXP(5)  *CMYP(5)  »IU(5) 
0X(52)  »DDX(54)  »DY(102) 

SNB ( 52)  * STB (52)  *UK(52*3) 

***  BLANK  COMMON 
DIMENSIONED  VARIABLES 
Z 

PK 

YY*  XX 

DDX*  DDY 

AMX*  AIX*  U» 

TAU*  JPM 

UL*  PL* 

XP*  YP*  CMXP* 


V (2502)  »P(2502)  » 

UPM(52)  » 

»  YAMC  ( 102) »  SIGCU02)* 


►  JK(5)  » 

►DDY ( 104)  » 

►  VK (52*3)  *RH0(52*3) 


variables 


XX 

DDY 

AIX* 

JPM 

PL* 

YP* 


non-oimensioned 


CMXP* 

variables 


CMYP 


COMMON 

AID 

» AMMV 

*ammy 

*AMPY 

»  Amur 

» AMUT  * AMVR  * 

1AMVT 

*DELEB  »DELER 

*DELET 

,delm 

*dtodx 

*DXYMIN*EAMMP  *EAMPY  • 

2E 

»ERDUMP* I 

*IWS 

*13  . 

*  J 

»K 

»KA  * 

3KB 

»LL  *MD 

»ME 

,MZT 

*nerr 

»NK 

*NPRINT* 

4NR 

* NRZ  vNULLE 

»PIDTS 

*siemin*snr 

*SNT 

*STR  * SOLID  » 

5SUM  *TESTRH»TWOPI 

6WFLAGL  *  WFLAGP 

*URR 

»WS  ' 

*  wsa 

*  WSB 

*WSC  » WFLAGF » 

COMMON  LAST 


***  THE  FOLLOWING  EQUIVALENCES  DEFINE  STORAGE  FOR 
X(0)»  Y(0) *  OX(O) *  DT(0) 

EQUIVALENCE  (XX(2)*  X(l))  » ( YY  <  2 ) *  Y(l))» 

1  (DDX (2) *DX( 1) )  * (0DY(2) »DY(1) ) 

***  SPECIAL  EQUIVALENCES  FOR  PH2  ONLY 

EQUIVALENCE  (UL*FLEFT>*  (UL( 103) * YAMC) * 

(PL*GAMC*PR>* 


(UL( 103) * YAMC) * 
(PL(103) *SXGC) 


***  SPECIAL  EQUIVALENCES  FOR  PH3  ONLY 
EQUIVALENCE  (UL»RSN)* 


(UL»RSN) * 
(PL»RST) * 

(P( 157) * VK) * 
(P(365)*STB>* 


( P » UK ) * 

(P(31 3) *SN8) * 
(P(417) »RHO) 


***  SPECIAL  EQUIVALENCES  FOR  EDIT 
EQUIVALENCE  (PR(1)*  IJ)»  (PR(6>*  JK>»  (UL(103) *CRAD) 
***  Z-STORAGE  EQUIVALENCES 
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o  o  o  o  o  o  o 


c 


EQUIVALENCE 

<Z( 

1) #PROB 

)»(Z( 

2) 

» cycle 

)# 

1<Z(  3)  »DT 

)»(Z( 

4) 

#  NUMSP 

)  # 

(Z  ( 

5) » 

NFRELP) 

»(Z( 

6) » 

NDUMP7) 

» 

2 (Z(  7)  # ICSTOP 

)  #  (Z( 

8) 

»PIDY 

)  # 

(Z( 

9) 

» TOPMU 

)  #  ( z 

10) 

»RTMU 

)  # 

3(Z(  U),STK1 

)»(Z( 

12) 

#NUMREZ) » 

(Z( 

13) 

»ETH 

)  #  (Z( 

14) 

» UN14 

)  # 

4(Z(  1S)#RHINIT 

)#(Z( 

16) 

#PROjI 

) » 

(Z( 

17) 

»KUNIT 

)  #  ( z 

18) 

9  XMAX 

)» 

5 ( Z (  19) »NZ 

)#(Z( 

20) 

»NREZ 

)  # 

(Z( 

21) 

» AMDM 

)#(Z( 

22) 

#  UVMAX 

)  # 

6 ( Z(  23) »UN23 

)»  (Z ( 

24) 

#DMIn 

)  # 

(Z( 

25) 

»JSTR 

)#(Z( 

26) 

»dtna 

)  # 

7(Z(  27)»CVIS 

)»(Z( 

28) 

#STK2 

)» 

(Z( 

29) 

»STEZ 

)»(Z( 

30) 

»NC 

)# 

8 ( Z (  31 ) »UN31 

)#(Z( 

32) 

» NRC 

) » 

( z 

33) 

»IMAX 

) » (Z( 

34) 

» IMAXA 

)# 

9(Z(  35)»JMAX 

)#(Z( 

36) 

» JMAXA 

)» 

<Z( 

37) 

»KMAX 

)»(Z( 

38) 

#KMAXA 

) 

EQUIVALENCE 
1(Z(  39) »BOTM 

)#(Z( 

40) 

»BOTmV 

)  # 

(Z( 

41) 

»NUMSPT 

)»(Z( 

42) 

»CZERO 

)  # 

2 ( Z (  43 ) » NUMSCA 

)#(Z( 

44) 

#PRLIM 

) » 

(Z( 

45) 

»PRDELT 

)»(Z( 

46) 

»PRFACT) 

EQUIVALENCE 
1(Z(  47) » 11 

)#(Z( 

48) 

#12 

) » 

<Z( 

49) 

» IPCYCL 

>*(Z( 

50) 

»TSTOP 

)  # 

2 ( Z (  51 ) rRHOFIL 

)#<Z( 

52) 

♦TARGV 

)  # 

(Z( 

53) 

»N3 

)#(Z( 

54) 

9 IVARDY) 9 

3(Z(  55) » VT 

)#(Z( 

56) 

» N6 

)» 

(Z( 

57) 

»rtm 

) » (Z( 

58) 

»rtmv 

)# 

4 ( Z (  59) #UN59 

)#(Z( 

60) 

#  N10 

) » 

(Z( 

61) 

#  Nil 

)#(Z( 

62) 

» gamma 

)  9 

5 ( Z(  63) »TOPM 

)r(Z( 

64) 

»BOTmU 

) » 

(Z( 

65) 

#SN 

)  9  (Z( 

66) 

»topmv 

)# 

6(Z(  67) »PRY80T 

)#(Z( 

68) 

»PRYTOP)» 

(Z( 

69) 

»PRXRT 

)»(Z( 

70) 

»CYCPH3) 9 

7(Z(  71)»REZFCT 

)»(Z( 

72) 

#  targi 

) » 

(Z( 

73) 

»PROJU 

)  9  (Z( 

74) 

rBBOUND) » 

8 ( Z (  75) »EVAP 

)#(Z( 

76) 

#ECK 

>* 

(Z  ( 

77) 

»NECYCL 

)  #  (Z( 

78) 

9  II 

)» 

9 ( Z <  79) » JJ 

)#  (Z( 

60) 

»NMP 

)» 

(Z( 

81) 

r  Y2 

)»(Z( 

82) 

»EZPH1 

) 

EQUIVALENCE 
1(Z(  63) » IVARDX 

)» <Z< 

84) 

»T 

>» 

(Z( 

85) 

»nmpmax 

)»(Z( 

86) 

»PMIN 

)9 

2 (Z(  67)# INTER 

)»(Z( 

68) 

»taybot>» 

<z< 

89) 

rTAYTOP 

)»(Z( 

90) 

»UN90 

)9 

3 ( Z (  91 ) » MC 

)»(Z( 

92) 

»MR 

)# 

(Z( 

93) 

»MZ 

)»(Z( 

94) 

#MB 

) 

EQUIVALENCE 

1(Z(  95) »REZ 

)»(Z( 

96) 

»N0DUMP) » 

(Z( 

97) 

»UN97 

)»(Z( 

98) 

»UN96 

)» 

2 ( Z(  99 ) » UN99 

) » (z(ioo) 

»EVAPM 

)  * 

(Z(101) 

»EVAPEN 

) 9 (Z ( 102) 

» EVAPMU) * 

3 ( Z ( 1 0  3 ) »EVAPMV 

) » (Z ( 104) 

»EZPh2 

) » 

<  Z ( 105) 

»SNL 

) # ( Z ( 106) 

#STL 

)» 

4 (Z (107 ) #  TAXRT 

) # (Z(1U8) 

»MSYmBL)» 

( 2  < 109) 

#UN109 

)#(Z(U0) 

»roeps 

)  9 

5 ( Z  < 111 ) »RHINI 

) # (Z( 112) 

» VINl 

)» 

( Z (113) 

»FINAL 

) 9  (Z(114) 

»FRSTD 

)» 

6 ( Z ( 1 15 ) » RHOZ 

) » ( Z ( 116) 

»ESA 

)» 

(Z ( 117) 

»esez 

) 9 (Z( 118) 

»esb 

)  9 

7(Z(119) »ESCAPA 

) # (Z( 120) 

»ESESP 

)» 

(Z(121) 

»ESESQ 

) 9  (Z( 122) 

»ESES 

)» 

6 ( Z ( 123 ) »ESALPH 

)» (Z(124) 

»ESB£TA)# 

( Z (125) 

» ESC APB 

) 9  (Z ( 126) 

»UN126 

)# 

9(Z(127)»SS1 

) # (Z ( 128) 

#SS2 

)» 

(Z ( 129) 

»UMIN 

) 9  (Z(130) 

» SS4 

) 

EQUIVALENCE 

1 (Z(131) tPRTIME 

) # (Z(132) 

»EOR 

) » 

(Z(133) 

#EOT 

) 9  ( Z ( 134 ) 

» EOB 

)  # 

2 (Z ( 135) »EMOR 

) # (Z( 136) 

»DXF 

)» 

(Z ( 137) 

»DYF 

) 9  (Z( 138) 

rRHOMIN) 9 

3 (Z(l39) # STAB) » 

(Z ( 140 ) # 

XIENrG) 

» 

(Z ( 141 ) 

» XKENRG 

)»  (Z(142)#XTENRG)» 

4(Z(143) »STT 

) # (Z ( 144) 

»DTMIN 

) » 

( Z ( 145) 

»trnsfc 

) » (Z( 146) 

»EMOT 

)# 

5(Z(147) » JPROJ 

)» (Z(148) 

»cnaut 

)* 

(Z ( 149) 

»bbar 

) » (Z(150) 

»EMOB 

) 

END  OF  COMMON 


DIMENSION  PROPI (50) »  AMK(l5)»  QK(15>»  TAB(15)»  CRAD(52) 

C  ***  SPECIAL  EQUIV.  FoR  EDIT 

EQUIVALENCE  (UL»PROPI)#  (PL<51)»  AMK) »  (PL(66)»QK)# 

1  <PL(81)»TAB> 

EQUIVALENCE  (PR(1)  »TIETAR) » (PR(2) rTKETAR)  »  (PRO)  »TETAR  )» 

1  (PR(4)»TARMAS)#(PR(5)»TARMV  ) » (PR(6) »TARMVP> » 


2  (PR (7) *RAMOMA) » (PR(8) »PRAMOA) *  <PR(9) *TIEPRO) * 

3  (PR(10)*TKEPRO)»(PR(11),TEPRO)*(PR(12)»PRMAS)* 

4  (PR ( 13) »PRMV  ),(PR(14),PRMVP)*(PR(15)*RAM0MB)* 

5  (PR ( 16) »PRAMOB) 

C 

C  . 

c 

C  ***  ENERGY  SUM  (ESUM)  AND  RELATIVE  ERROR  IN  SUM  (RELERR) 

C  COMPUTED.  ECK  IS  LARGEST  ERROR  COMPUTED  AND  ON  PRINT 

C  CYCLES  IS  PRINTED  AND  COMPARED  TO  DMIN*  MAXIMUM 

C  ALLOWABLE  ERROR. 

C 

C  . 

ESUM=0. 

00  1U  K=2»KMAX 

1 0  ESUM=ESUM+ AMX (K ) * ( • 5* (U (K ) **2+V ( K ) **2 ) ♦AIX (K) ) 

RELEKR= ( ESUM-ETH ) /ETH 
IF  (ABS(RELERR) .LT.AtiS(ECK) )  GO  TO  20 
ECKSRELERR 
NECYCL=NC 
20  CONTINUE 

C  ***  NERR  =  1  WHEN  ERROR  CALLS  EDIT. 

IF  (NERR.EQ.l)  GO  TO  150 

C  . . .  •  •  • 

c 

C  ***  NPRINT  =  1  WHEN  EDIT  IS  CALLED  TO  DO  AN  INTERMEDIATE 

C  PRINT.  SKIP  TESTS  ON  TIME  TO  STOP*  PRINT*  REZONE* ETC. 

C  WHICH  ALREAOY  HAVE  BEEN  DONE  FOR  THIS  CYCLE. 

C 

C  . . .  • 

IF  (NPRINT.EQ.1)  GO  TO  190 
C  ***  13=1  SIGNALS  A  SHORT  PRINT 

13=1 

C  ***  IF  THIS  IS  FIRST  CYCLE  OF  RUN*  WFLAGF=1. 

IF  (WFLAGF.GT.O.)  GO  TO  120 

C  ***  IS  THIS  THE  TIME  OR  CYCLE  TO  STOP  EXECUTION 

IF  ( ICSTOP.LE.NC.AND. ICSTOP.GT.O)  GO  TO  140 
IF  (T*(1.+ROEPS).GE.TSTOP.AnD.TSTOP,GT.O.)  GO  TO  140 
C  ***  SHOULD  THE  GRID  BE  REZONED 

IF  ( (REZ.NE.0. .AND.REZFCT .NE.O, . AND.NUMREZ.GT.O) .OR.SS4.NE.O. )  GO 
1T0  145 
C 

ASSIGN  414  TO  LOCA 
ASSIGN  110  TO  LOCB 

C  . .  . . . . . 

c 

C  ***  ARE  WE  PRINTING  ON  TIME  OR  CYCLE  INTERVALS 

C 

C  . . . .  •  •  • 

40  IF  (PRDELT.NE.O.)  GO  TO  50 
IF  ( IPCYCL.NE.O)  GO  TO  100 
GO  TO  420 

C  ***  PRINTING  ON  TIME.  IS  IT  TIME  TO  PRINT 

50  IF  (TM1.+ROEPS)  .GE.PRTIME)  GO  TO  70 

C  ***  NO.  BUT  WILL  NEXT  CYCLE  BYPASS  THE  PRINT  TIME 

IF  (PRTIME.GE.T+DT)  GO  TO  60 
DT=PRTIME-T 
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t 


*1 


< 

DTNA=DT 

60  60  TO  LOCA#  (414.412) 

•C  ***  YES#  IT  IS  TIME  TO  PRINT.  nPRINT=1  FLAGS  THIS  AS  A 

C  PRINT  CYCLE. 

70  NPRlNT-1 

C  ***  AVOID  TRUNCATION 

T=PRTIME 

C  ***  IS  IT  TIME  TO  RESCALE  PRINT  INTERVAL 

IF  (T*(l.+ROEPS) .LT.PRLIM.OR.NUMSCA.LE.O)  GO  TO  80 
C  ***  CHANGE  PRINT  INTERVAL  AND  THE  TIME  FOR  THE  NEXT 

C  RESCALING. 

PRDELT=PRDELT*PRFACT 

PRLIM=PRLIM*PRFACT 

NUMSCA=NUMSCA-1 

C  *"*  DEFINE  TIME  FOR  NEXT  PRINT. 

80  PRTIME=T+PRDELT 

I WS= ( PRT IME+ . 5*PRDELT ) /PRDElT 
WS=IwS 

PRTIME=WS*PRDELT 

C  ***  WILL  WE  BYPASS  TIME  To  PRInT 

IF  (PRTIME.GE.T+DT)  60  TO  90 
C  ***  YES#  ADJUST  DT 

ut=pktime-t 

dtna=ot 

90  GO  TO  LOCB#  (110#412) 

C  ***  PRINTING  ON  CYCLES.  IS  THIS  A  PRINT  CYCLE 

100  IF  (MOD(NC#IPCYCL).NEVO>  GO  TO  LOCA#  (414# 412) 

C  ***  YES.  NPRINT  =  1  FLAGS  THIS  AS  A  PRINT  CYCLE. 

NPRINT=1 

c  ***  is  this  The  cycle  to  rescale  print  interval 

IF  (NC. LT.PRLIM.OR.NUMSCA.LE.O)  GO  TO  LOCB#  (110*414) 

C  ***  YES.  MULTIPLY  NUMBER  OF  CYCLES  BETWEEN  PRINTS  BY  PRFACT 

IPCYCL=INT(PRFACT)*IPCYCL 
PRLIM=PRFACT*PRLIM 
NUMSCA=NUMSCA-1 
GO  TO  LOCB.  (110.412) 

C  ***  TEST  FOR  SHORT  OR  LONG  PRINT 

C  ***  NUMSP  COUNTS  NUMBER  OF  SHORT  PRINTS  SINCE  LAST  LONG 

C  PRINT.  NUMSPT  COUNTS  NUMBER  OF  PRINTS  SINCE  LAST 

C  TAPE  DUMP. 

110  NUMSP=NUMSP+1 

NUMSPT=NUMSPT+1 

IF  ( NUMSP. NE.NFRELP)  GO  TO  130 
NUMSP=0 

C  ***  13=11  SIGNALS  A  LONG  PRINT 

120  13=11 

C  ***  PRINT  OF  RESTART  CYCLE  WILL  BE  SHORT  IF  PK(3).LT.«1. 

IF  (PK(3).LT.-1..AND.WFLAGF.GT,0.)  13=1 
130  IF(NUMSPT.NE.NDUMP7)  GO  TO  190 
GO  TO  150 

C  ***  SET  WFLAGL=1.  TO  SAY  THIS  IS  LAST  CYCLE  OF.  RUN 

140  WFLAGL=1. 

145  13=11 

NPRINT=1 

NUMSP=0 


C 

C 
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c  ***  tape  dump 

c 

c  . 

150  NUMSPT=0 

IF  (NFRELP.EQ.NDUMP7)  NUMSRrO 
IF  (NODUMP. NE.O)  GO  TO  170 
BACKSPACE  KUNIT 
IF (FKSTD.GT.O. )  GO  TO  155 
BACKSPACE  KUNIT 
155  WS-555.0 

WRIT*.  (KUNIT)  WS»CYCLE»N3 
WRITE  (KUNIT)  (Z (L) *L=1*MZT) 

WRITE  (KUNIT)  (U(K) *V(K) »AMx (K) »AIX(K) *P(K) »K=1*KMAXA> 

WRITE  (KUNIT)  X(0)*(X(K)*TAU(K)»JPM(K>*K=1*IMAX) 

WRITE  (KUNIT)  (Y(K) »K=0»JMAX> 

C  ***  ARE  TRACER  POINTS  BEING  GENERATED 

IF  (Y2.GT • (“1. ) )  GO  TO  160 

WRITE  (KUNIT)  ((XP(I#J)#YP(I#J)#I=1»II)#J=1,JJ) 

160  WRITE  (KUNIT)  (DX ( I ) » 1=1 * IMaX) 

WRITE  (KUNIT)  (DY(I> »I=1*JMaX) 

WS=666.0 

WRITE  (KUNIT)  WStWS'WS 
FRSTO  =0. 

END  FILE  KUNIT 
170  CONTINUE 

C  ***  ERDUMP=1.  WHEN  ERROR  CALLS  EDIT  FOR  A  TAPE  DUMP  ONLY 

IF  (ERDUMP.GT.O.)  RETURN 

c  . 

c 

C  ***  COMPUTE  AND  PRINT  ENERGY*  MASS  AND  MOMENTUM  TOTALS. 

C 

C  . . 

C  ***  INITIALIZE  PR  ARRAY*  TEMPORARY  STORAGE  FOR  ENERGY* MASS 

C  AND  MOMENTUM  TOTALS  PRINTED  OUT. 

190  DO  200  1=1*16 

200  PR(I)=0. 

C 

C  RAf10MA=RADTAL  MOMENTUM  ABOVE  UPROJ 

C  RAMOMA=RADTAL  MOMENTUM  BELOw  UPROJ 

C  PRAMOA=POSITIVE  RADIAL  MOMENTUM  ABOVE  JPROJ 

C  PRAMOB=POSITIVE  RADIAL  MOMENTUM  BELOW  JPROJ 

C 

IF  (JPROJ. NE.O)  GO  TO  205 
N=2 

GO  TO  220 

205  N=IMAX*JPROJ+l 

DO  210  K=2*N 
WS=AMX(K> 

PRMAS=PRMAS+WS 

TIEPRO=TILPRO+WS*AIX(K) 

TKEPRO=TKEPRO+  .  5«WS*  ( U  ( K > ♦♦2+V  ( K )  **2 ) 

WSA=WS*V(K) 

prmv=prmv+wsa 

IF  (WSA.GT.O.)  PRMVP=PRMVP*wSA 
R AMOMBSRAMOMB+AMX ( K ) *U ( K ) 

IF  (U(K).GT.O.)  PRAMOB=PRAMOB+AMX(K)*U(K) 

210  CONTINUE 
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N=N+i 

220  00  250  K=N#  KMAX 

hS=AMX(K) 

tarmas=takmas+ws 

TIETAR=TIETAR+WS*AIX(K) 

TKETAR=TKETAR+.5*WS*(U(K)**2+V(K)**2) 

wSA=wS*V(K) 

TARMV=TARMV+WSA 

IF  (wSA.GT.O.)  TARMVP=TARMVP+WSA 
R AMOMA=RAMOMA+AMX (K) *U(K) 

IF  (U(K).GT.O.)  PRAMOA=PRAMOA+AMX(K)*U(K) 

250  CONTINUE 

TETAR=TIETAR+TK£TAR 

tepro=t:epro+tkepro 

00  240  J=l»8 
PR(J+16)=PR(J)+PR(J+8) 

240  CONTINUE  „ 

IF  (IMAX.GT.I)  GO  TO  260 
C 

c  ***  if  doing  a  i-d  problem  divide  totals  by  nz  where 

C  NZ=4** (NUMBER  OF  TIMES  THE  GRID  HAS  BEEN  REZONED* ) 

C 

PROPI (1)=ETH/NZ 
PROPI (2)=ECK/NZ 
PR0P1(4)=EZPH1/NZ 
PR0Pl(5)=EZPH2/NZ 
PROPI ( 6) -BBOUND/NZ 
DO  250  J-l »24 
250  PROPI (J+6)=PR(J)/NZ 
PROPI ( 31 )=BOTM/NZ 
PROPI (32) =RTM/NZ 
PROPI (33)=T0PM/NZ 
PROPI (34) =EVAPM/NZ 
PROPI (35) =EMOB/NZ 
PROPI ( 36) =EMOR/NZ 
PROPI (37)=EM0T/NZ 
PROPI ( 38 ) -E VAPEN/NZ 
PROPI (39)-B0TMU/NZ 
PROPI ( 40 )=RTMU/NZ 
PROPI (41)=T0PMU/NZ 
PROPI (42)=EVAPMU/NZ 
PROPI (43)=B0TMV/NZ 
PROPI (44) =RTMV/NZ 
PROPI ( 45) =TOPMV/NZ 
PROPI (46) =EVAPMV/NZ 
PROPI ( 47) =EOB/NZ 
PROPI  (48UE0R/NZ 
PROPI ( 49 )=EOT/NZ 

WRITE  (6# 520)  PR0B#T»NC#PR0PI(1)»PR0PI(2)#NECYCL# (PROPI ( J) * J=4^6) 
WRITE  (6» 530)  (PROPI ( J) * Jz7*49) 

GO  TO  270 

260  WRITE  (6»520)  PR0B»T»NC»ETH»ECK*NECYCL»EZPH1»EZPH2»BB0UND 

WRITE  (6#530)  ( (PR (U) » J=l»24) >BOTM#rTM#TOPM»EVAPM»EMOB»EMOR*EMOT>E 
1VAPEN»BOTMU»RTMU»TOPMU#EVAPmU»BOTMV»RTMV»TOPMV»EVAPMV»EOB»EOR*EOT) 
270  WRITE  (6»570)  ( JPM( I ) #  1=1* II ) 

C  ***  ENERGY  TOTALS  STORED  FOR  LATER  USE  IN  TRACER  POINT 

C  PLOTS. 
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XIENKG=PR<17) 

XKENRG=PR(18) 

XTENRG=PR(19> 

C  . 

NKT  =  12 

C  ***  TABS  ARE  TANGENT  ALPHAS 


TABU) 

s  0.02 

TABU) 

=  0.04 

TAB (3) 

=  0.06 

TAB (4) 

=  0.08 

TAB (5) 

=  0.10 

TABU) 

=  0.15 

TABU) 

=  0.20 

TABU) 

=  0.25 

TAB (9) 

=  0.30 

TAB ( 10 ) 

=  0.40 

TAB (11) 

=  0.50 

TAB( 12) 

=  1.00 

c 

NK1  =  NKT+2 

DO  275 

1-1 » NKl 

AMK ( I ) 

=  0. 

PK  ( I ) 

=  0. 

27b 

GiK(I) 

=  0. 

C 

DO  280 

K52.KMAXA 

IF(AMX(K) )440»28Q»276 

276 

I=NK1 

IF(V(K) ) 279 #279.277 

277 

wSA  s  ABS(U(K))/V(K) 

C 

C  ***  SEARCH  FOR  TAN  AnGLE  MADE  BY  VELOCITY  VECTOR  OF  CELL. 

C 

DO  278  I=1»NKT 
IF(TABd)-WSA)  278»279»  279 

278  CONTINUE 
I=NK+1 

C  ***  SUM  MASS  BETWEEN  ANGLES. 

279  AMK(I)  =  AMK(I)  ♦  AMX(K) 

C  ***  SUM  RAOIAL  MOMENTA  BETWEEN  ANGLES. 

PK ( I )  =  PK ( I )  ♦  U(K) *AMX(K) 

C  ***  SUM  AXIAL  MOMENTA  BETWEEN  ANGLES. 

QK ( I )  S  QK ( I )  +  V (K ) *AMX ( K ) 

280  CONTINUE 
WRITE<6»605) 

WRITE (6. 610) ( AMK ( I ) » 1  =  1 » NKl ) 

WRITE (6r 615) (PK(I) #  I=1»NK1) 

WRITE (6 » 620) (QK(I) »  I=1»NK1) 

IF  (NUMSPT.EQ.O)  WRITE (6» 540)  NC 
C 

C  ***  ARE  TRACER  POINTS  BEING  GENERATED 

IF  (Y2.GT • (•!.) )  GO  TO  305 


C 

C  ***  PRINT  TRACER  POINT  COORDINATES  IN  CM. 

C 

C  . . . 
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i; 

i, 


:| 


c 

C  ***  IJ,  JK  =  THE  I  and  J  OF  THE  CELL  THE  TRACER  POINT 

C  ORIGINATED  IN  .  (TRACER  POINTS  CHANGE  POSITION  IN 

C  XP  AND  YP  ARRAYS  WHEN  THEY  ARE  WEEDED  OUT 

C  DURING  REZONE.) 

C 

I J  ( N )  =2**  <  NRZ+1 )  *  <  1-1  M-l 
JK(N)=2**(NRZ*l)*(J-i)+l 
IF  (N.LT.5)  GO  TO  300 

WRITE  (6*500)  (IJ(M),JK(M)»CMXP(M)*CMYP(M)»M=1,N) 

N=0 

300  CONTINUE 

IF  (N.EQ.O)  GO  TO  305 

WRITE  (6*500)  (IJ(M) »JK(M) »CMXP(M) »CMYP(M) »M=1»N) 

305  IF  ( IMAX.EQ. 1 )  GO  TO  360 

C  . 

c 

C  PRINT  SYMBOLIC  CONTOUR  MAPS  OF  COMPRESSION*  PRESSURE* 

C  VELOCITY,  AND  INTERNAL  ENERGY  UNLESS  DOING  A  1-D 

C  PROBLEM. 

C 

C  . . . 

C 

CALL  MAP 

C 

C  . . . 

c 

C  *♦*  COMPUTE  CRATER  DEPTH  AND  VOLUME.  AID  SUMS  DEPTH. 

C 

C  . . . . . . 

WRITE(6»490) 

AID  =  0. 

C  ♦**  START  AT  AXIS 

DO  330  I  =1*11 
CRAD(I)  =  .5*DX(I)+X(I-1) 

PL( I )  =0. 

UL ( I )  =  0. 

DO  320  J  =1*12 
K=(J-1)*IMAX  +1+1 
C  ***  WS  IS  COMPRESSION 

WS  =  AMX(K)/(TAU(I)*DY(J)*RHOZ) 

XFCWS.LT. (.99))  GO  TO  310 


WRITE  (6*580) 

N-Q 

D 0  3U0  J=1»JJ 
DO  3U0  I=1»II 

IF  (XP(I*J).LE.0..AND.YP(1*J).LE.0.)  GO  TO  300 
IP=INT(XP(I*J) ) 

JP=INT (YP( I , J) ) 

KK=JP*IMAX+IP+2 

IF  ( AMX(KK) ,GT .0. )  GO  70  290 

XP(I,J)=0. 

YP( I » J) -0. 

GO  TO  300 
290  N=N+1 

CMXP(N)=X(IP)4-DX(IP+1)*(XP(I,J)-INT(XP(I»J)) ) 
CMYP(N)=Y(JP)4-DY(JP+1)*(YP(I»J)-INT(YP(I»J))) 
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GO  TO  325 

310  AID  =  AID  ♦  l.-WS 

C  ***  NOT  AT  BOTTOM  OF  CRATER  YET 

320  CONTINUE 
325  IAIO  =  I NT (AID) 

C  ***  UL(I)  IS  CM.  DEPTH  OF  CRATER  IN  COLUMN  I 

C  ***  PLd)  IS  CELL  DEPTH  OF  CRATER  IN  COLUMN  I 

UL( I)  =  Y(IAID)  ♦  DY<IAID+l)*<AlD-FLOAT<IAID>)  -  Y(UPROJ) 
IF(UL(I).GT*O..OR.UL(I).LT.Q.)  PL<I)  =  AID 
AID  =  0. 

330  CONTINUE 

C  ***  PRINT  CRATER  DEPTHS 

JO  340  1=1*11 

IF(UL(I) *LT*0**OK*UL(I) .GT.O.)  60  TO  335 
GO  TO  340 

335  WRITE(6»495)  I*  PL(I)»  CRAD(I)»  UL<I) 

340  CONTINUE 

C  ***  COMPUTE  CRATER  VOLUME  AND  VOLUME  OF  HEMISPHERE  WITH 

C  f^ADlUS-UL(  1)  • 

WSB-0. 

DO  345  I=1»U 
IF < UL < I > .LT»0.)  60  TO  350 
C  ***  WSB  6IVES  CRATER  VOLUME 

wSB  =  UL( I ) *TAU( I )  4- WSB 
345  CONTINUE 
350  CONTINUE 

C  ***  PRINT  CRATER  VOLUME  ONLY  WHEN  GREATER  THAN  ZERO 

IF (WSb.GT.O. )  60  TO  355 
GO  TO  360 

C  ***  WSC  GIVES  VOLUME  OF  HEMISPHERE 

355  *SC52.0944*(UL(1))**3 

WRITt<6'498)  WSBf  WSC 

C  . 

c 

C  ***  SHORT  PRINT  MEANS  13=1  AND  PROPERTIES  ARE  PRINTED  ONLY 

C  FOR  CELLS  IN  FIRST  COLUMN*  LONG  PRINT  MEANS  13=11  AND 

C  PROPERTIES  ARE  PRINTED  FOR  ALL  CELLS  IN  ACTIVE  GRID. 

C 

C  . * . 

360  DO  410  151*13 
KSPACE=G 
WFLAGPsl. 

J=I2+1 

K=l2*IMAX+l+I 
DO  400  L=i*I2 
J=J-1 
K=K-IMAX 

365  IF  ( AMX(K) )  440*390*370 

370  IF  (WFLAGP.EQ.O*)  GO  TO  380 

WRITE  (6*550)  I*X<I)*DX(I) 

WFLAGP=0. 

380  WS=AMX(K)/(TAU(I)*OY(J}) 

WSA=WS/RHOZ 
V»SC=P  ( K ) 

WRITE  (6*510)  J*U(K)*V(K)*WSC'AMX(K)*WS*AIX(K)*WSA*Y(J) 

KSPACE=0 
GO  TO  400 
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390 

KSPACE=KSPACE+1 

IF  (KSPACE.GT.l)  GO  TO  400 
WRITE  (6»560) 

400 

CONTINUE 

410 

r 

CONTINUE 

# 

W 

IF  (NPRINT.EG.l)  GO  TO  412 
ASSIGN  412  TO  LOCA 

ASSIGN  412  TO  LOCB 

IF  (PRDELT.NE.O.)  GO  TO  50 
GO  TO  100 

c 

***  CHECK  ON  SIZE  OF 

ENERGY  DISCREPANCY 

412 

IF  (ABS(ECK) .GT.DMIN)  GO  TO 

430 

C 

***  IF  LAST  CYCLE#  REWIND  TAPE 

414 

IF  ( wFLAGL.EG.O. )  GO  TO  416 
REWIND  7 

416 

wFLAGP=0. 

WFLAGF=0. 

c  ***  nerr=i  when  edit  is  called  by  error. 

IF  (NERR.EQ.l)  return 


C  ***  SHOULD  GRID  BE  ReZONED  ON  THIS  CYCLE 

IF  ( (REZ.NE.O • .AND.REZFCT . NE»0 • • AND.NUMREZ.GT.O) . OR.SS4.NE.O. )  GO 
1TO  419 

C  . .  .  . 

c 

RETURN 

c 

c  . . 

c 

C  ***  R  E  Z  0  N  E 

C 

c  . . . . 

419  CALL  REZONE 

C  ***  MUST  CALL  CDT  TO  RECALCULATE  PRESSURES 

TNOW=T 
DTNOw=DT 
REZ=0. 

SS4=Q. 

CALL  CDT 
T=TNOW 
DT=DTNOW 
DTNASDT 

NUMREZ=NUMREZ-1 

C 

C  ***  NREZ  =  NUMBER  OF  REZONES  ALLOWED  (INPUT  VALUE  OF  NUMREZ) 

C  NUMREZ  =  NUMBER  OF  REZONES  ALLOWED  MINUS  THE  NUMBER 

C  OF  REZONES  PERFORMED  SINCE  T=0. 

C 

NRZ=NREZ-NUMREZ 

C  ***  NZ  USED  IN  PRINTOUT  OF  TOTALS  FOR  1-D  PROBLEMS 

NZ=4.**NRZ 
C 

GO  TO  145 

C  . . . . . 

c 

C  ***  ERROR  CONDITIONS 

C 
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***  PRINT  DELTA  NOT  SPECIFIED  IN  INPUT 


f 

c 

••  c 


f*20 

NK=4U 

GO  TO  450 

C 

*** 

ENERGY  CHECK 

430 

NKS412 

GO  TO  450 

C 

***  NEGATIVE  MASS 

440 

NK=365 

450 

NR=5 

CALL  ERROR 

c 

C  FORMATS 

C 

490  FORMAT  ( 1H0 #  17X  # 35HDEPTH  OF  CRATER  MEASURED  FROM  JPR0J//12X»1HI»5X 
1»18HJ  OF  CRATER  BOTTOM# 12X# 1HR# 11X» 17H0EPTH  IN  CM.  D(I)//) 

495  FORMAT  ( 113# 9X# 0PF6. 1 » 13X » IpElO . 4 » 9X » 1PE10.4) 

496  FORMAT  ( //6X # 13HCRATER  VOLUME# 11X»43HCRATER  VOLUME  BASED  ON  (2/3) 
1*  PI  *  D(1)**3/7X#1PE10.4#26X»1PE10.4) 

500  FORMAT  (5(14# 14# 1P2E9.2) ) 

510  FORMAT  ( 14 # IX# 1P2E14.6# 3E15.6#E14. 6#E15.6»E14*6) 

52  0  FORMAT ( 8H1PR0BLEM  #  6X » 4HT I ME  #  8X  #  5HC YCLE » 3X » 1 3HT0T . EN . THEOR . 3X » 

1  19HMAX • REL • ERROR-CYCLE #  3X » 18H I E  SET  TO  ZER0-PH1»3X» 

2  18HIE  SET  TO  ZER0-PH2#3X» 12HPLASTIC-W0RK/1F8.4»2X» 1PE13.7# 

3  3X»I4»4X»1PE13.7»3X»1PE13.7»1X#I4»6X#1PE13.7#8X#1PE13.7»6X# 

4  1PE13.7/) 

530  FORMAT  (18X»2HIE» 14X#2HKE#7x»l3HT0T.EN.  (SUM) »7X»4HMASS» 12X»2HMV»8 
1X»12HMV (POSITIVE) #8X#2HMU»8X>12HMU(P0SITIVE)/11H  J.GT. JPROJ# 1P8E15 

2.7/11H  J. LE . JPROJ #  1P6E15. 7/14X » 12H - - - -  12H - —-•»—-# 

33X » 1 2H— — - — —  #  3X » 12H - — — - -# 3X » 1 2H - - - - » 3X » 12H - 

4- - - »3X » 12H— - - - 3X  *  12H - - #3X/7H  TOTALS#  4X » IP 

58E15 . 7///9H  BOUNDARY » 9X » 6HB0TT0M » 9X » 5HRIGHT » 10X » 3HT0P #  8X » 12R4EV APO 
6RATED4//9H  MASS  0UT»2X»1P4E15.7/11H  ENERGY  OUT# 1P4E15.7/7H  m£)  OUT# 
74X»1P4E15.7/7H  MV  0UT»4X» 1P4E15. 7//11H  WORK  DONE  »1P3E15.7//) 

540  FORMAT  (1H0//21H  TAPE  7  DUMP  ON  CYCLEI5////) 

550  FORMAT  (1H  ///4H  I  =13# 6X#6hR( I )  =F12.3»6X#7HDR ( I )  =E14.7//3H  J8X 
1 # 1HU13X# 1HV13X# 3H  P  12X# 3HAMX12X# 3HRH0UX»3HAIX12X#4HC0MP11X»2H  Z/ 
2) 

560  FORMAT  (1H0) 

570  FORMAT  (//22H  J  OF  PRESSURE-MAXlMUM/(25I5) ) 

580  FORMAT (//103H  TRACER  POINTS  -  INITIAL,  LOCATION  IN  CELL  COORDINATES 
1  (I#J)  -  CURRENT  LOCATION  IN  CM.  COORDINATES  (X# Y)//  5(4H  I#3X» 

21HJ»5X» 1HX»8X» 1HY » 3X) ) 

605  FORMAT (//41H  ANGULAR  DISTRIBUTION  OF  MASS  AND  MOMENTA/130H  TAN  U/V 
1  0-.Q2  .02-. 04  .04-.06  .06-. 08  .08-#10  .10-.15  .15-. 20  .20- 

2.25  .25-. 30  .30-. 40  .40-. 50  .50-1.0  1.0-UP  V.LE.O./) 

610  FORMAT (5H  MASS# 14(1X» 1PE8.2) ) 

615  FORMAT (5H  MU  »14(1X»1PE8.2) ) 

620  FORMAT (5H  MV  #14(1X#1PE8.2) ) 

ENO 


1 


c 

c 


c 

c 


c 


c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 


SUBROUTINE.  MAP 


DIMENSION 

AMX (2502) 

*  AlX (2502) *U(2502)  • 

V(2502) 

#P(2502)  » 

1 

X  ( 52) 

# XX (54)  »TAU(52)  » 

JPM (52) 

t 

2 

Y ( 102) 

» YY (104)  »FLEFT ( 102) 

»  YAMC (102) »  SIGC (102) 

3 

GAMC (102) 

t 

4 

PK ( 15) » 

Z  <  150 )  » 

5 

XP(26»  51) 

»YP(26*5l)» 

6 

PL(204) 

# UL (204)  »PR (204)  # 

7 

RSN (52) » 

RST(52)» 

8  v. 

CMXP(5) 

#CMYP(5)  *IJ(5)  » 

JK  ( 5) 

» 

9 

DX(52) 

t DDX ( 54)  »DY ( 102)  • 

DDY (104) 

» 

S 

SNB ( 52 ) 

» STB ( 52)  »UK (52#  3)  » 

VK(52#3) 

•RH0(52»3) 

***  B  L  A  N 

K  COMMON 

DIMENSIONED  VARIABLES 

COMMON 

Z 

COMMON 

PK 

COMMON 

YY# 

XX 

COMMON 

DDX  t 

DDY 

COMMON 

AMX# 

AIX#  U* 

V# 

P 

COMMON 

TAU# 

JPM 

COMMON 

UL» 

PL# 

COMMON 

XP# 

YP»  CMXP* 

CMYP 

non-dimensioned  variables 

COMMON 

aid 

* AMMV  »AMMY  »AMPY 

» AMUR 

»amut  »amvr 

1AMVT  ,DELEB  »DELER  »DELET  ,DELM  »DTODX  >DXYMIN»EAMMP  »EAMPY  * 
2E  » EROUMP* I  » IMS  »I3  »J  »K  »KA  * 


3KB  »LL  »MD  »ME  »MZT  »NERR  »NK 

4NR  »NRZ  »NULLE  *PIOTS  #SIEMIN»SNR  »SNT 

5SUM  »TESTRH»TWOPI  »URR  »WS  »WSA  »WSB 
6wFLAGL»  WFLAGP 
COMMON  LAST 


*NPRINT» 

#STR  »SOLID  » 
*WSC  »WFLA6F  * 


***  THE  FOLLOWING  EQUIVALENCES  DEFINE  STORAGE  FOR 
X(0)»  Y ( 0 ) »  DX(0)»  DY(O) 

EQUIVALENCE  (XX (2) »  X ( 1 ) )  » (YY(2) »  Y( 1) ) » 

1  (DDX(2) »DX( 1) )  »(DDY(2> »OY(l)) 


***  SPECIAL 
EQUIVALENCE 

1 

***  SPECIAL 


EQUIVALENCES  FOR  PH2  ONLY 

(UL*FLEFT) t  (UL( 103) rYAMC) * 

(PL»GAMC»PR> »  (PL(103) tSIGC) 

EQUIVALENCES  FOR  PH3  ONLY 


EQUIVALENCE 

1 

2 

3 


(UL»RSN) » 
(PLrRST) t 
(P( 157) » VK) » 
(P(365) ,STB) t 


(PrUK) » 
(P(3l3) »SNB>  * 
(P(417)»RH0) 


***  SPECIAL  EQUIVALENCES  FOR  EOIT 


EQUIVALENCE  (PR(1)»  IJ)»  (PR<6)»  JK> t  (ULC103) »CRAD) 


***  Z-STORAGE  EQUIVALENCES 
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o  o  o  o  o  o  o 


EQUIVALENCE 
1(Z(  3)  »DT 

2 (Z(  7)  # ICSTQP 

3  ( Z  (  11)#STK1 

4 <Z<  15) #RHINIT 
5(Z(  19) » NZ 
6 ( Z (  23 ) » UN23 
7<Z<  27)#CVIS 
8 ( Z (  31 ) #UN31 
9(Z(  35) » JMAX 
EQUIVALENCE 
1(Z(  39) »BOTM 
2(Z(  43) »NUMSCA 
EQUIVALENCE 
1 ( Z (  47 ) » I 1 
2 ( Z (  51 ) »RHOFIL 
3(Z(  55) » VT 
4 ( Z (  59) #UN59 
5(Z (  53) » TOPM 
6(Z(  67) »PKYBOT 
7(Z(  71) ?KEZFCT 
8 (Z(  75)  e EVAP 
9(Z(  79) » JJ 
EQUIVALENCE 
1(Z(  83) » IVARDX 
2 (Z (  87) » INTER 
3(Z(  91 ) »MC 
EQUIVALENCE 
1(Z(  95) »REZ 
2 (Z(  99) »UN99 
3(Z(103) »EVAPMV 

4  ( 1  ( 107) #TAXRT 
EiZ(lll) »RHINI 
6 (Z ( 115) »RHOZ 
7(Z(119)»ESCAPA 
8 ( Z ( 123 ) » ESALPH 
9 ( Z ( 127) #SS1 

EQUIVALENCE 
1 (Z(131 ) »PRTIME 
2 (Z ( 135) »EM0R 
3(Z(139) »STAQ) » 
4(Z(143) »STT 
5(Z(147) » JPROU 


(Z( 

1) 

»prob  ')»(Z( 

2) 

»CYCLE 

)  # 

(Z( 

4) 

#  NUMSP 

)*(Z( 

5)  » 

NFRELP) » (Z( 

6) » 

SIDUMP7) 

» 

(Z( 

8) 

#PIDY 

)  # 

(Z< 

9) 

» TOPMU  )»(Z( 

10) 

»RTMU 

)  # 

(Z( 

12) 

»  numrez) » 

(Z( 

13) 

#ETH  )»(Z( 

14) 

» UN14 

) » 

(Z( 

16) 

#  PROjI 

)» 

(Z( 

17) 

#KUNIT  )#(Z( 

18) 

» XMAX 

)  » 

(Z  ( 

20) 

#NREZ 

) » 

(Z( 

21) 

» AMDM  )»(Z( 

22) 

#  UVMAX 

)  # 

(Z( 

24) 

f  DMIn 

)» 

(Z( 

25) 

» JSTR  )»(Z( 

26) 

»DTNA 

) » 

(Z( 

28) 

»STK2 

)» 

(Z( 

29) 

»STEZ  )»(Z( 

30) 

»NC 

) » 

(Z( 

32) 

»NRC 

)» 

(Z( 

33) 

» IMAX  )»(Z( 

34) 

» IMAXA 

) » 

(Z( 

36) 

» JMAxA 

)» 

(Z( 

37) 

»KMAX  )»(Z( 

38) 

»KMAXA 

) 

(Z( 

40) 

♦BOTmV 

)» 

(Z( 

41) 

#  NUMSPT) » (Z ( 

42) 

#CZERO 

)  » 

(Z( 

44) 

»PRLIM 

)» 

(Z( 

45) 

#PRDELT)»(Z( 

46) 

»PRFACT 

) 

(Z( 

48) 

#12 

)» 

(Z( 

49) 

» IPCYCL) # (Z( 

50) 

#TSTOP 

) » 

(Z( 

52) 

#  TAR&V 

)» 

<Z( 

53) 

»N3  )»(Z( 

54) 

» IVARDY 

) » 

(Z( 

56) 

»N6 

)» 

(Z( 

57) 

»RTM  )#(Z( 

58) 

#RTMV 

)  » 

(Z( 

60) 

#  N10 

)# 

(Z( 

61) 

#  Nil  )#(Z( 

62) 

» GAMf^A 

)  » 

(Z( 

64) 

»BOTmU 

)# 

(Z( 

65) 

»SN  )»(Z( 

66) 

» TOPJijV 

) » 

(Z( 

68) 

#PRYTOP)» 

(Z( 

69) 

»PRXRT  )»(Z( 

70) 

»CYCf»H3 

) » 

(Z( 

72) 

# TARSI 

)» 

(Z « 

73) 

» PROUU  )#(Z( 

74) 

#BBOUNO 

)  » 

(Z( 

76) 

»ECK 

)» 

(Z( 

77) 

#NECYCL)»(Z( 

78) 

»II 

) » 

(Z( 

80) 

#NMP  ' 

) » 

(Z( 

81) 

»Y2  )#(Z( 

82) 

»EZPH1 

) 

(Z( 

84) 

#  T 

)» 

(Z( 

85) 

» NMPMAX) » (Z( 

86) 

»PMIN 

)  # 

(Z( 

88) 

» TAYbOT) » 

(Z( 

89) 

» TAYTOP) » (Z( 

90) 

#UN90 

)  # 

(Z( 

92) 

#  MR 

)# 

(Z  ( 

93) 

#MZ  )#(Z( 

94) 

#MB 

) 

(Z( 

96) 

» NODUMP)# 

(Z  ( 

97) 

#UN97  )»(Z( 

98) 

»UN98 

)  # 

(Z(100) 

» EVApM 

)» 

(Z(101) 

#  EVAPEN) » (Z ( 102) 

rEVAPMl 

1)  # 

(Z ( 104) 

»EZPh2 

)# 

(Z ( 105) 

* SNL  ) »  (Z ( 106) 

»STL 

)  # 

(Z( 108) 

»MSYmBL)» 

(Z ( 109) 

» UN109  )#(Z(110) 

»ROEPS 

)  # 

( Z  ( 1 12 ) 

»VINI 

)» 

<Z( 113) 

#FINAL  ) » (Z( 114) 

»FRSTO 

)  # 

(Z(116) 

» ESA 

)# 

(Z  < 117) 

#ES£Z  ) » (Z( 118) 

»ESB 

)  # 

( Z( 120 ) 

#ESE5P 

)» 

( Z ( 121 ) 

» ESESQ  )  #  (Z( 122) 

»ESES 

)  # 

(Z ( 124) 

#ESBETA)» 

(Z ( 125) 

»ESCAPB) # (Z(126) 

#UN12G 

)# 

( Z ( 128) 

#SS2 

)» 

( Z ( 129) 

# UMIN  ) »  (Z ( 130 ) 

#SS4 

) 

(Z( 132) 

»EOR 

)# 

(Z ( 133) 

»EOT  )»(Z(134) 

»EOB 

)# 

(Z ( 136) 

r  DXF 

)# 

(Z(137) 

#DYF  )  # (Z(138) 

»RHOMIN) » 

Z ( 140 ) # 

XIENRG)# 

( Z ( 141 ) 

# XKENR6) #  (Z(142)»XTENRG) 

(Z( 144) 

» DTMIN 

)# 

(Z( 145) 

#TRNSFC)»(Z(146) 

»EMOT 

)  # 

(Z ( 148) 

#cnaut 

)» 

( Z ( 149) 

#BBAR  ) » (Z( 150) 

»EMOB 

) 

END  OF  COMMON 


DIMENSION  PROP (52) #  WSMAX(5)»  VALUE(41> 

C  ***  SPECIAL  EQUIVALENCE  FOR  MAP 

EQUIVALENCE  (UL#PROP)»  (Ul(52) rWSMAX) #  (UL(157) » VALUE  ' 

C 

DIMENSION  ALE (41) 

DATA  ALE/  2H  .»2H  -#2H  A»2H  B»2H  C»2H  D#2H  E#2H  F# 


22 


■c  w  t\>  ^  oj  w  h 


C 

C 

c 

c 

1 


2 

C 

3 


4 

C 


6 

C 


8 

C 


10 

C 

c 


2h 

G#2H 

H»2H 

I  #2H 

J»2H 

K  »2h 

L»2H 

M,2H 

N»2H 

O' 

2H 

P  #  2H 

Q»2H 

R»2H 

S»2H 

T»2H 

U»2H 

V»2H 

W » 2H 

X' 

i  2H 

Y » 2H 

Z»2H 

♦  #2H 

* »  2h 

1 » 2H 

2»2H 

3»2H 

4 » 2H 

5' 

2H 

b#2H 

7»2H 

8»2H 

9 » 2H 

0  »2H 

/ 

DIMENSION  XUM ( 41 ) 

uATA  XUM/ 

2H 

• » 2H 

-#  2H-A  »  2H-B » 2H-C » 2H-D » 2H-E  *  2H-F » 

2H-G » 2H-H » 2H-1 » 2H- J » 2  H"K » 2H-L » 2H-M » 2H-N » 2H-0 » 

2H-P » 2H-Q , 2H-R » 2H-S » 2  H"T »  2H-U » 2H-V » 2H-W » 2H-X » 

2H-Y  »  2H-2 » 2H-+ » 2H-* » 2H-1 » 2H-2 » 2H-3 » 2H-4  #  2H-5» 
2H-b»2H-7»2H-8#2H-9»2H-0»2H  / 

IDL=MIN0( Il»b4> 

JDL=12 

IF  (NC.NE.O)  GO  TO  1 
IDL=MIN0 ( I MAX » 54) 

JDL-JMAX 

***  FIND  MAXIMUM  VALUE  IN  ACTIVE  GRID  OF  EACH  PROPERTY 

***  COMPRESSION 
*SMIN=10E20 
WSMAX ( 1 ) =0 • 

DO  2  J=1»JDL 

DO  2  I=1»IDL 

K- ( J-l ) + I MAX* I+1 

IF(AMX(K) .EQ.O.)  GO  TO  2 

COMP  =  AMX(K)/(DY(J)*TAU(I)*RHOZ) 

WSMAX ( 1 )  =  AMAXKWSMAX(l)  »COMP> 

WSMIN  =  AMIN1  (WSMIN'COMP) 

CONTINUE 

IF (WSMAX (1) #GT»WSMlN)  GO  TO  3 
wSMlN  -  0  • 

***  PRESSURE 
wSMAX(2>=0, 

DO  4  J=l»JOL 
DO  4  I=1»IDL 
K=(  J-1)*XMAX+I+1 

WSMAX ( 2 )  =  AMAX1(WSMAX(2)»ASS(P(K))> 

***  RADIAL  VELOCITY 
WSMAX ( 3) =0. 

DO  6  J=l»JOL 
DO  6  I=l»IOL 
K=(J-1)*IMAX+I+1 

rtSMAX (3)  =  AMAX1 ( WSMAX ( 3) »  AbS (  U  (  K )  )  ) 

***  AXIAL  VELOCITY 
WSMAX<4)=0, 

DO  8  J=1»JDL 
DO  8  1  =  1 » IDL 
K=(J-1)*IMAX+I«-1 

WSMAX ( 4 )  =  AMAX1(WSMAX(4) »ABS( V<K) ) ) 

***  SPECIFIC  INTERNAL  ENERGY 
WSMAX ( 5) =0 • 

DO  10  J=l»JOL 
DO  10  1  =  1  * IDL 
K=( J-l ) *IMAX+I+1 

WSMAX (5)  =  AMAX1(WSMAX(5) »AB$< AlX(K) ) ) 

***  STORE  INFORMATION  TO  BE  PLOTTED  IN  PROP  ARRAY 


C  A  ROW  AT  A  TIME. 

C 

NPROP  =  i 

C  ***  COMPRESSION 

J=JDL 

MS=MSYMBL+1 

WRITE(6»500) 

15  DO  20  I=1»IDL 

K=(J-1)*IMAX+I+1 

20  PROP ( 1 )  =  AMX(K)/(TAU(I)*DY(J)*RHOZ) 

GO  TO  110 
C 

C  ***  PRESSURE 

C 

30  J=JDL 

MS=MSYMBL 

WRITE(6»510) 

35  DO  40  1=1* IDL 

K=(  J<*1 )  *IMAX+I  +  1 
40  PROP ( 1 )  =  P(K) 

GO  TO  110 
C 

C  ***  RADIAL  VELOCITY 

C 

50  J=JDL 

WRITE (6*  520 ) 

55  DO  60  1=1' IDL 

K={ J-l ) ♦IMAX+I+1 
60  PROP ( I )  =  U(K) 

GO  TO  110 
C 

C  ***  AXIAL  VELOCITY 

C 

70  J=JDL 

WRITE (6» 530) 

75  DO  80  1=1' IDL 

K=(J-1)*IMAX+I+1 
80  PROP ( I )  =  V (K) 

GO  TO  110 
C 

C  ***  SPECIFIC  INTERNAL  ENERGY 

C 

90  J=JDL 

WRITE(6»540) 

95  DO  100  1=1' IDL 

K=(J-1)*IMAX+I41 
100  PROP(I)  =  AIX(K) 

C 

C  ***  WHEN  PRINTING  FIRST  (TOP)  ROW  OF  MAP»  COMPUTE 

C  SCALE  FACTOR  AND  PRINT  KEY. 

110  IF(J.LT.JDL)  GO  TO  300 
C 

IF(WSMAX (NPROP) .GT.O. )  GO  TO  180 
C  ***  SKIP  CALCALATION  OF  SCALE  FACTOR 

GO  TO  300 
C 

C  ***  COMPUTE  SCALE  FACTOR  AND  PRINT  MAXIMUM  VALUE  OF 

2k 


C  EACH  SYMBOL  USED 

C 

JdO  SCALE  =  WSMAX(NPROP)/FLOAT(MS) 

IF  ( (AINT(SCALE*1000. ) ) .LT. (SCALE*1000.) )  GO  TO  190 
00  TO  200 

190  SCALE  =  AINT(SCALE*1000.+1>/1000. 

200  CONTINUE 
C 

IF(NPROP.EQ.I)  GO  TO  220 
VALUE ( 1 )  =  0. 

VALUE ( 2 )  £  SCALE/10. 

DO  210  I=1#MS 

210  VALUE ( 1+2)  =  FLOAT(I)*SCALE 

GO  TO  240 

C  ***  VALUES  FOR  COMPRESSION  MAP 

220  VALUE ( 1 )  =  WSMIN 
DO  2,50  1  =  1 » MS 

230  VALUE ( I  +  l )  =  FLOAT ( I ) *SCALE 
C  ***  PRINT  DEFINITIONS  OF  MAP  SYMBOLS 

240  ILIM1  =  1 
ILIM2  =  10 
MSP=MSYMBL  +  2 

250  IF  ( MSP.LT . ILIM2)  ILIM2  =  MSP 
IF  (NPROP.NE.l)  GO  TO  260 
WRITE <6» 550)  (ALE(I) #I=ILIM1»ILIM2) 

WRITE (6» 560)  (VALUE(I) »I=ILIM1»ILIM2) 

GO  TO  270 

260  WRITE ( 6 »  570 )  < ALE ( I ) » I=ILIMl » ILIM2) 

WRITE(6.580)  (VALUE(I) »I=ILIM1»ILIM2) 

270  IF  (MSP.EQ.ILIM2)  GO  TO  280 
ILIM1=ILIM2+1 
ILIM2=ILIM2>10 
GO  TO  250 
280  WRITE (6»  590) 

C 

C  ***  ASSIGN  APPROPRIATE  SYMBOL  TO  EACH  CELL  IN  ROW  J. 

C 

300  DO  370  I=1»IDL 

K=( J-1)*IMAX+I+1 
IF  ( AMX(K) *GT *0. )  GO  TO  310 
MA  =  41 
GO  TO  360 

310  IF ( NPROP.EQ. 1 )  GO  TO  340 

IF(AoS(PROP(I) ) .GT.O. )  GO  To  320 
MA  =  1 
GO  TO  360 

320  IF ( At)S ( PROP (I)).GT.VALUE(2)  )  GO  TO  330 

MA  =  2 
GO  TO  360 

330  FLOTMA  =  ABS(PROP( I ) ) /SCALE  +  2. 

MA  =  INT (FLOTMA) 

IF(FLOTMA.GT.AlNT(FLOTMA) )  ma=ma+i 
MA  =  MAXO (MA»3) 

GO  TO  360 

C  ***  DEFINE  MA  FOR  COMPRESSION  MAP 

340  IF(PR0P(I) .GT. WSMIN)  GO  TO  350 
MA=1 
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GO  TO  360 

350  FLOTMA  =  ASS (PROP(I)) /SCALE  ♦  1. 

MA  =  IMT (FLOTMA) 

IFCFLOTMA.GT.AINT(FLOTMA) )  MA  =  MA+1 
MA  =  MAXO ( MA » 2) 

C  ***  STORE  CHARACTER  TO  BE  PLOTTED  FOR  CELL  K 

360  PR  ( I )  =  ALE(MA) 

IF(PROP(I) .LT.O.)  PR  ( I )  =  XUM(MA) 

C  ***  END  OF  I-LOOP 

370  CONTINUE 

C  ***  PRINT  J  ROW  OF  MAP 

IF(M0D(J»5) .NE.O)  GO  TO  380 
WRIT£(6»600)  J,  (PR ( I )  * I— 1 * IDL ) 

GO  TO  390 

380  WRITE (6*610 )  (PR(I)»  I=1»ID|_) 

390  J=J-1 

C  ***  HAVE  WE  REACED  BOTTOM  ROW 

IF(J.EQ.O)  GO  TO  395 
GO  TO  (15#35»55#75»95) »NPROP 
C  ***  PRINT  AND  LABEL  X-AXIS  OF  MAP 

395  PR ( 1 )  =  ALE (29) 

wRITt (6*600)  J*  (PR(1) *1=1* IDL) 

WRITE(6*620)  (I*  1=0 » IDL*  5) 

C 

NPROP  =  NPROP  ♦  1 
GO  TO  (400 *30 *50 *70 *90 *400) » NPROP 
C 

400  RETURN 

C  ***  FORMATS 

500  FORMAT ( 1H1 » 4X » 15HC0MPRESSI0N  //) 

510  FORMAT (1H1»4X». 5HPRESSURE  //) 

520  FORMAT ( 1H1 » 4X» 15HRADIAL  VELOCITY//) 

530  FORMAT (1H1*4X* 15H AXIAL  VELOCITY  //) 

540  F ORMaT ( 1H1 * 4X * 24HSPECIFIC  INTERNAL  ENERGY//) 

550  FORMAT ( 16H  SYMBOL  *l0(4X.A2,4X) ) 

560  FORMAT ( 16H  MAXIMUM  VALUE  » 10 (F6.3»4X) > 

570  FORMAT ( 16H  SYMBOL  » 10(3X» A2*5X) ) 

580  FORMAT (16H  MAXIMUM  VALUE  *1P10E10.2) 

590  FORMATS//) 

600  FORMAT ( 110* 2H  I»54A2) 

610  F0RMAT(10X*2H  I»54A2) 

620  FORMAT(U2*10I10////) 

END 
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REQUIREMENTS  FOR  ENLARGING  THE  GRID 


The  RIM  code,  as  it  is  listed  in  this  report,  can  calculate  at  most 
2500  cells.  The  number  of  rows  (JMAX)  cannot  exceed  100,  and  the  number 
of  columns  (IMAX)  cannot  exceed  50.  To  increase  the  grid  size  the  user 
needs  only  to  redimension  most  of  the  arrays  in  Blank  Common  and  redefine 
some  of  the  equivalences.  Given  IMAX  and  JMAX,  the  parameter  definitions 
belo^:  show  how  the  dimensions  and  equivalences  should  be  redefined. 

PARAMETERS;  IJD  =  IMAX  *  JMAX  +  2 
ID  =  IMAX  +  2 
IDP  =  IMAX  +  4 
JD  =  JMAX  +  2 
JDP  =  JMAX  +  4 
ITP  =  (IMAX  +  2)/2 
JTP  =  (JMAX  +  2  )/2 
JD2  =  2  *  (JMAX  +  2) 

ID3  =  3  *  IDP 
ID4  =  6  *  IDP 
ID5  =  9  *  IDP 
ID6  =  12  *  IDP 

DIMENSIONS; 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

# 


AMX(IJD),  AIX(IJD),  U(IJD),  V(UT>),  P(UT>), 
X(lD),  XX(lDP),  TAU(ID),  JIM(lD), 

Y(JD),  YY(JDP),  FLEFT(JD),  YAMC(JD),  SIGC(JD), 
GAMC(JD), 

PK(15),  Z(150), 

XP(ITP,  JTP),  YP(lTP,JTP), 

PL(JD2),  UL( JD2) ,  PR(JD2), 

RSN(ID),  RST(ID), 

CMXP(5),  CMYP ( 5 ) >  IJ(5),  JK(5), 

DX(ID),  DDX(IDP),  DY(JD),  DDY(JDP), 

snb(id),  STB(ID),  UK(ID,3),  VK(ID,3),  RH0(ID,3) 
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EQUIVALENCES: 


(UL,FLECT),  (UL(JDP),YAMC) 
(PL,GAMC,PR),  (PL(JDP),SIGC) 

(P(ID3),VK),  (P(ID4),SI1B), 
(P(ID5),STB),  (P(n)6),RHD) 

(PR(1),IJ),  (PR(6),JK),  (UL(JDP),CRAD) 
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